Analyses

Set-up

Packages

library(brms)
library(ggplot2)
library(lme4)
library(performance)
library(see)
library(sjmisc)
library(tidyverse)

options(
  digits = 3
)
set.seed(170819)

Data

d <- read_csv("data/data.csv")
Rows: 960 Columns: 464
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr   (45): persistence, anonymity, topic, QUESTNNR_T1, MODE_T1, DE03_09_T1,...
dbl  (407): case, n_OE, n_Words, mean_Words, number_messages, CASE_T1, REF_T...
lgl    (8): SERIAL_T1, IV01_01_T1, K001_CP_T1, K001_T1, MAILSENT_T1, SERIAL_...
dttm   (4): STARTED_T1, LASTDATA_T1, STARTED_T2, LASTDATA_T2

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# same as above; but original file wording
# d <- read_csv("data/DataAggregated_T1T2_costsbenefits.csv")
# load("data/image.RData")

# recode to make as sum coding
d$anonymity_dev <- factor(d$anonymity)
contrasts(d$anonymity_dev) <- contr.sum(2)
d$persistence_dev <- factor(d$persistence)
contrasts(d$persistence_dev) <- contr.sum(2)

d <- d %>% 
  rename(
    group = roles,
    op_expr = n_OE
  )

Descriptives

Let’s inspect distribution of opinion expressions.

ggplot(d, aes(op_expr)) +
  geom_histogram(binwidth = 1)

Looks like a zero-inflated poisson distribution. Confirms our preregistered approach to analyze data using zero-inflated Poisson approach.

nrow(d[d$op_expr == 0, ]) / nrow(d)
[1] 0.214

Overall, 21% of participants without any opinion expressions.

Let’s look at distribution of experimental groups.

d %>% 
  select(persistence, anonymity) %>% 
  table
           anonymity
persistence high low
       high  240 240
       low   240 240

Distribution among groups perfect.

d %>% 
  select(topic) %>% 
  table
topic
  climate    gender migration 
      320       320       320 

Distribution of topics also perfect.

Let’s check if groups are nested in topics:

is_nested(d$topic, d$group)
'f2' is nested within 'f1'
[1] TRUE

Indeed the case.

We first look at the experimental group’s descriptives

d %>% 
  group_by(persistence) %>% 
  summarize(op_expr_m = mean(op_expr))
# A tibble: 2 × 2
  persistence op_expr_m
  <chr>           <dbl>
1 high             9.26
2 low              9.29

Looking at persistence, we see there’s virtually no difference among groups.

d %>% 
  group_by(anonymity) %>% 
  summarize(op_expr_m = mean(op_expr))
# A tibble: 2 × 2
  anonymity op_expr_m
  <chr>         <dbl>
1 high           9.03
2 low            9.51

People who with less anonymity communicated more. But the difference isn’t particularly large.

d %>% 
  group_by(persistence, anonymity) %>% 
  summarize(op_expr_m = mean(op_expr))
`summarise()` has grouped output by 'persistence'. You can override using the
`.groups` argument.
# A tibble: 4 × 3
# Groups:   persistence [2]
  persistence anonymity op_expr_m
  <chr>       <chr>         <dbl>
1 high        high           9.27
2 high        low            9.25
3 low         high           8.79
4 low         low            9.78

Looking at both groups combined, we see that low anonymity and low persistence created highest participation. But differences among groups aren’t large.

d %>% 
  group_by(group) %>% 
  summarize(
    anonymity = anonymity[1],
    persistence = persistence[1],
    topic = topic[1],
    op_expr_m = mean(op_expr)
    ) %>% 
  rmarkdown::paged_table()

Looking at the various individual groups, we do see some difference. Generally, this shows that individually, communication varied. But not so much systematically across experimental groups.

d %>% 
  group_by(topic) %>% 
  summarize(op_expr_m = mean(op_expr))
# A tibble: 3 × 2
  topic     op_expr_m
  <chr>         <dbl>
1 climate        9.63
2 gender         9.96
3 migration      8.22

Looking at topics specifically, we also see that there’s some variance.

Manipulation Check

Let’s see if respondents indeed perceived the experimental manipulations.

model_pers <- lm(
  per_persistence ~ persistence_dev,
  d
)

summary(model_pers)

Call:
lm(formula = per_persistence ~ persistence_dev, data = d)

Residuals:
   Min     1Q Median     3Q    Max 
-3.194 -0.777 -0.110  0.806  3.223 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        2.9854     0.0378    79.1   <2e-16 ***
persistence_dev1  -1.2085     0.0378   -32.0   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1 on 703 degrees of freedom
  (255 observations deleted due to missingness)
Multiple R-squared:  0.593, Adjusted R-squared:  0.592 
F-statistic: 1.02e+03 on 1 and 703 DF,  p-value: <2e-16
# report::report_text(model_pers)

The experimental manipulation worked.

model_anon <- lm(
  per_anonymity ~ anonymity_dev,
  d
)
summary(model_anon)

Call:
lm(formula = per_anonymity ~ anonymity_dev, data = d)

Residuals:
   Min     1Q Median     3Q    Max 
-3.425 -0.216 -0.216  0.575  3.784 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      2.8204     0.0348    81.0   <2e-16 ***
anonymity_dev1  -1.6042     0.0348   -46.1   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.924 on 703 degrees of freedom
  (255 observations deleted due to missingness)
Multiple R-squared:  0.751, Adjusted R-squared:  0.751 
F-statistic: 2.12e+03 on 1 and 703 DF,  p-value: <2e-16
# report::report_text(model_anon)

The experimental manipulation worked.

Bayesian mixed effects modeling

We analyze the data using Bayesian modelling.

We use deviation/sum contrast coding (-.1, .1). Meaning, contrasts measure main effects of independent variables.

Fixed effects

We preregistered to analyze fixed effects. However, we did not explicate if we were to model zero inflation separately. Below hence two options.

fit_fe_1 <- 
  brm(
    op_expr ~ 
      1 + persistence_dev * anonymity_dev +
      (1 | topic/group)
    , data = d
    , chains = 4
    , cores = 4
    , iter = 8000
    , warmup = 2000
    , family = zero_inflated_poisson()
    , control = list(
      adapt_delta = .95
      , max_treedepth = 12
      )
    , save_pars = save_pars(all = TRUE)
    )
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
using SDK: ‘MacOSX14.4.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Warning: There were 767 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

Only some convergence issues. Let’s inspect model.

plot(fit_fe_1, ask = FALSE)

Looks quite good!

Let’s look at results.

summary(fit_fe_1)
Warning: There were 767 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
 Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev + (1 | topic/group) 
   Data: d (Number of observations: 960) 
  Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
         total post-warmup draws = 24000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.26      0.28     0.01     1.05 1.00     1615      753

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.40      0.05     0.32     0.51 1.00     2590     6421

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           2.38      0.19     1.93     2.78 1.00
persistence_dev1                   -0.00      0.06    -0.12     0.13 1.01
anonymity_dev1                      0.00      0.06    -0.11     0.12 1.00
persistence_dev1:anonymity_dev1     0.03      0.06    -0.08     0.15 1.00
                                Bulk_ESS Tail_ESS
Intercept                           1318      344
persistence_dev1                    1207      369
anonymity_dev1                      3229     6017
persistence_dev1:anonymity_dev1     3081     6134

Further Distributional Parameters:
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi     0.21      0.01     0.19     0.24 1.00     6898     6452

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

No significant effect emerged.

Let’s inspect ICC

var_ratio_fe <- performance::variance_decomposition(
  fit_fe_1
  , by_group = TRUE)
var_ratio_fe
# Random Effect Variances and ICC

Conditioned on: all random effects

## Variance Ratio (comparable to ICC)
Ratio: 0.43  CI 95%: [-0.14 0.73]

## Variances of Posterior Predicted Distribution
Conditioned on fixed effects: 29.50  CI 95%: [13.80 58.84]
Conditioned on rand. effects: 51.60  CI 95%: [46.86 56.86]

## Difference in Variances
Difference: 22.06  CI 95%: [-7.11 38.06]

42.8177583 percent of variance in opinion expressions explained by both topics and groups.

Let’s look at exponentiated results for interpretation.

broom.mixed::tidy(fit_fe_1) |> 
  mutate(
    estimate_exp = exp(estimate),
    conf_low_exp = exp(conf.low),
    conf_high_exp = exp(conf.high)
    ) %>% 
  mutate(
    across(
      c(
        estimate
        , conf.low
        , conf.high
        , estimate_exp
        , conf_low_exp
        , conf_high_exp
        ),
      ~ round(.x, 2)
    )
  ) %>% 
  select(
    -effect, -component, -group, -std.error
  ) %>% 
  rmarkdown::paged_table()
Warning in tidy.brmsfit(fit_fe_1): some parameter names contain underscores:
term naming may be unreliable!

Again, of course no significant effect. Main effects appear trivial. Interaction effect somewhat larger. Let’s visualize results to see what this exactly means.

p <- plot(
  conditional_effects(
    fit_fe_1
    ), 
  ask = FALSE,
  plot = FALSE
  )

p_anon <- p[["anonymity_dev"]]
p_int <- p[[3]]

p_pers <- 
  p[["persistence_dev"]] +
  xlab("Persistence") +
  ylab("Opinion expression") +
  scale_x_discrete(
    breaks = c(-0.5, .5),
    labels = c("Low", "High")
    ) +
  scale_y_continuous(
    limits = c(5, 15)
  )

p_anon <- 
  p[["anonymity_dev"]] +
  xlab("Anonymity") +
  ylab("Opinion expression") +
  scale_x_discrete(
    breaks = c(-0.5, .5),
    labels = c("Low", "High")
    ) +
  theme(
    axis.title.y = element_blank()
    ) +
  scale_y_continuous(
    limits = c(5, 15)
  )

p_int <- 
  p[[3]] +
  xlab("Persistence") +
  ylab("Opinion expression") +
  scale_x_discrete(
    breaks = c(-0.5, .5),
    labels = c("Low", "High")
    ) +
  scale_color_discrete(
    labels = c("Low", "High")
    ) +
  guides(
    fill = "none",
    color = guide_legend(
      title = "Anonymity"
      )
    ) +
  theme(
    axis.title.y = element_blank()
    ) +
  scale_y_continuous(
    limits = c(5, 15)
  )

plot <- cowplot::plot_grid(
  p_pers, p_anon, p_int, 
  labels = c('A', 'B', "C"), 
  nrow = 1,
  rel_widths = c(2, 2, 3)
  )

plot

ggsave("figures/results.png", plot, width = 8, height = 4)

Shows that there are no main effects. There seems to be a (nonsignificant) interaction effect. In high persistence environment, anonymity is conducive to communication; in low it’s the opposite.

Let’s look at posteriors

p_1 <- 
  pp_check(fit_fe_1) + 
  labs(title = "Zero-inflated poisson")
Using 10 posterior draws for ppc type 'dens_overlay' by default.
p_1

The actual distribution cannot be precisely reproduced, but it’s also not too far off.

Random effects

We said we’d explore random effects. Following the “keep it maximal” principle, it’d actually make sense to consider these analyses the best ones, as they can model how the experimental conditions affect the outcomes differentially depending on topic.

Let’s now model random effects.

fit_re_1 <- 
  brm(
    op_expr ~ 
      1 + persistence_dev * anonymity_dev +
      (1 + persistence_dev * anonymity_dev | topic) + 
      (1 | topic:group)
    , data = d
    , chains = 4
    , cores = 4
    , iter = 8000
    , warmup = 2000
    , family = zero_inflated_poisson()
    , control = list(
      adapt_delta = .95
      , max_treedepth = 15
      )
    , save_pars = save_pars(all = TRUE)
    )
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
using SDK: ‘MacOSX14.4.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Warning: There were 1281 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems

Shows some problems, but not all that many.

Let’s inspect model.

plot(fit_re_1, ask = FALSE)

Traceplots look alright.

Let’s look at results.

summary(fit_re_1)
Warning: There were 1281 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
 Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev + (1 + persistence_dev * anonymity_dev | topic) + (1 | topic:group) 
   Data: d (Number of observations: 960) 
  Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
         total post-warmup draws = 24000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
                                                      Estimate Est.Error
sd(Intercept)                                             0.36      0.46
sd(persistence_dev1)                                      0.25      0.35
sd(anonymity_dev1)                                        0.24      0.34
sd(persistence_dev1:anonymity_dev1)                       0.41      0.46
cor(Intercept,persistence_dev1)                          -0.00      0.46
cor(Intercept,anonymity_dev1)                             0.01      0.46
cor(persistence_dev1,anonymity_dev1)                     -0.01      0.47
cor(Intercept,persistence_dev1:anonymity_dev1)            0.08      0.46
cor(persistence_dev1,persistence_dev1:anonymity_dev1)     0.02      0.46
cor(anonymity_dev1,persistence_dev1:anonymity_dev1)      -0.01      0.46
                                                      l-95% CI u-95% CI Rhat
sd(Intercept)                                             0.01     1.60 1.00
sd(persistence_dev1)                                      0.01     1.21 1.00
sd(anonymity_dev1)                                        0.00     1.20 1.00
sd(persistence_dev1:anonymity_dev1)                       0.02     1.69 1.00
cor(Intercept,persistence_dev1)                          -0.83     0.83 1.00
cor(Intercept,anonymity_dev1)                            -0.83     0.83 1.00
cor(persistence_dev1,anonymity_dev1)                     -0.83     0.83 1.00
cor(Intercept,persistence_dev1:anonymity_dev1)           -0.78     0.86 1.00
cor(persistence_dev1,persistence_dev1:anonymity_dev1)    -0.83     0.84 1.00
cor(anonymity_dev1,persistence_dev1:anonymity_dev1)      -0.83     0.83 1.00
                                                      Bulk_ESS Tail_ESS
sd(Intercept)                                             6620    11752
sd(persistence_dev1)                                      7764    10211
sd(anonymity_dev1)                                        7280     9370
sd(persistence_dev1:anonymity_dev1)                       6024     5585
cor(Intercept,persistence_dev1)                          21881    16269
cor(Intercept,anonymity_dev1)                            22853    15065
cor(persistence_dev1,anonymity_dev1)                     19146    13364
cor(Intercept,persistence_dev1:anonymity_dev1)           22045    15452
cor(persistence_dev1,persistence_dev1:anonymity_dev1)    19117    15472
cor(anonymity_dev1,persistence_dev1:anonymity_dev1)      15719    17072

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.40      0.05     0.32     0.51 1.00     7776    12613

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           2.38      0.27     1.79     2.95 1.00
persistence_dev1                   -0.00      0.21    -0.41     0.44 1.00
anonymity_dev1                      0.00      0.22    -0.44     0.44 1.00
persistence_dev1:anonymity_dev1     0.02      0.33    -0.67     0.71 1.00
                                Bulk_ESS Tail_ESS
Intercept                          10960     9291
persistence_dev1                   10694     8366
anonymity_dev1                      9075     6054
persistence_dev1:anonymity_dev1     8583     5583

Further Distributional Parameters:
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi     0.21      0.01     0.19     0.24 1.00    32167    15762

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Again, virtually no main effect. Interaction effect larger, but also not even closely significant. Let’s see if the random effects model fits better

fit_fe_1 <- add_criterion(fit_fe_1, "loo", moment_match = TRUE)
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Warning: Found 13 observations with a pareto_k > 0.7 in model 'fit_fe_1'. With
this many problematic observations, it may be more appropriate to use 'kfold'
with argument 'K = 10' to perform 10-fold cross-validation rather than LOO.
fit_re_1 <- add_criterion(fit_re_1, "loo", moment_match = TRUE)
Error in validate_ll(log_ratios) : 
  All input values must be finite or -Inf.
Error in eval(expr, envir, enclos) : 
  Exception: model6c7925fe16d3__namespace::write_array: Cor_1 is not positive definite. (in 'anon_model', line 188, column 2 to column 66)
Error in eval(expr, envir, enclos) : 
  Exception: model6c7925fe16d3__namespace::write_array: Cor_1 is not positive definite. (in 'anon_model', line 188, column 2 to column 66)
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Warning: Found 17 observations with a pareto_k > 0.7 in model 'fit_re_1'. With
this many problematic observations, it may be more appropriate to use 'kfold'
with argument 'K = 10' to perform 10-fold cross-validation rather than LOO.
comp_1 <- loo_compare(fit_fe_1, fit_re_1, criterion = "loo")
comp_1
         elpd_diff se_diff
fit_fe_1  0.0       0.0   
fit_re_1 -1.0       2.4   

Although model comparisons showed that the model with random effects fitted better, the difference was not significant (Δ ELPD = -0.9779785, 95% CI [-5.7077202, 3.7517632]. Hence, for reasons of parsimony the model with fixed effects is preferred.

Hurdle

Let’s now estimate a fixed effects model with hurdles.

fit_hrdl_1 <- 
brm(
  bf(
    op_expr ~ 
      1 + persistence_dev * anonymity_dev +
      (1 | topic) + 
      (1 | topic:group),
    zi ~ 
      1 + persistence_dev * anonymity_dev +
      (1 | topic) + 
      (1 | topic:group)
  )
  , data = d
  , chains = 4
  , cores = 4
  , iter = 8000
  , warmup = 2000
  , family = zero_inflated_poisson()
  , control = list(
    adapt_delta = .95
    , max_treedepth = 15
    )
  )
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
using SDK: ‘MacOSX14.4.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Warning: There were 601 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems

Estimation works quite well.

Let’s inspect model.

plot(fit_hrdl_1, ask = FALSE)

Trace-plots look alright.

summary(fit_hrdl_1)
Warning: There were 601 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
 Family: zero_inflated_poisson 
  Links: mu = log; zi = logit 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev + (1 | topic) + (1 | topic:group) 
         zi ~ 1 + persistence_dev * anonymity_dev + (1 | topic) + (1 | topic:group)
   Data: d (Number of observations: 960) 
  Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
         total post-warmup draws = 24000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)        0.27      0.32     0.01     1.17 1.00     2245     1183
sd(zi_Intercept)     0.28      0.44     0.01     1.50 1.00     6240     5389

~topic:group (Number of levels: 48) 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)        0.40      0.05     0.32     0.50 1.00     5272    10476
sd(zi_Intercept)     0.24      0.14     0.01     0.53 1.00     5361     8349

Regression Coefficients:
                                   Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                              2.39      0.20     1.96     2.83 1.00
zi_Intercept                          -1.30      0.28    -1.74    -0.67 1.00
persistence_dev1                      -0.01      0.06    -0.12     0.11 1.00
anonymity_dev1                         0.00      0.06    -0.12     0.12 1.00
persistence_dev1:anonymity_dev1        0.03      0.06    -0.09     0.15 1.00
zi_persistence_dev1                    0.03      0.09    -0.15     0.20 1.00
zi_anonymity_dev1                      0.01      0.09    -0.17     0.18 1.00
zi_persistence_dev1:anonymity_dev1     0.01      0.09    -0.16     0.18 1.00
                                   Bulk_ESS Tail_ESS
Intercept                              2994     1307
zi_Intercept                           7245     4630
persistence_dev1                       4984     9340
anonymity_dev1                         5664     9023
persistence_dev1:anonymity_dev1        5982     9817
zi_persistence_dev1                   19443     9630
zi_anonymity_dev1                     22174    16988
zi_persistence_dev1:anonymity_dev1    10375     7268

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Same results, no main effects, slightly larger but still nonsignificant interaction effect.

Let’s inspect coefficients individually.

coefficients(fit_hrdl_1)
$topic
, , Intercept

          Estimate  Est.Error     Q2.5    Q97.5
climate   2.441330 0.09198736 2.268439 2.630276
gender    2.416212 0.09056169 2.244328 2.603659
migration 2.306579 0.09761436 2.106201 2.485353

, , zi_Intercept

           Estimate Est.Error      Q2.5     Q97.5
climate   -1.312747 0.1255348 -1.561256 -1.056106
gender    -1.331712 0.1254764 -1.587955 -1.087517
migration -1.340756 0.1259943 -1.602138 -1.098537

, , persistence_dev1

              Estimate  Est.Error       Q2.5     Q97.5
climate   -0.005103491 0.05928069 -0.1215837 0.1109697
gender    -0.005103491 0.05928069 -0.1215837 0.1109697
migration -0.005103491 0.05928069 -0.1215837 0.1109697

, , anonymity_dev1

              Estimate  Est.Error       Q2.5     Q97.5
climate   0.0009798597 0.05871981 -0.1159714 0.1169367
gender    0.0009798597 0.05871981 -0.1159714 0.1169367
migration 0.0009798597 0.05871981 -0.1159714 0.1169367

, , persistence_dev1:anonymity_dev1

            Estimate  Est.Error        Q2.5     Q97.5
climate   0.02947201 0.05854357 -0.08580015 0.1467483
gender    0.02947201 0.05854357 -0.08580015 0.1467483
migration 0.02947201 0.05854357 -0.08580015 0.1467483

, , zi_persistence_dev1

            Estimate  Est.Error      Q2.5     Q97.5
climate   0.02959441 0.08905165 -0.148826 0.2034417
gender    0.02959441 0.08905165 -0.148826 0.2034417
migration 0.02959441 0.08905165 -0.148826 0.2034417

, , zi_anonymity_dev1

             Estimate  Est.Error       Q2.5     Q97.5
climate   0.006613642 0.08968717 -0.1673214 0.1817552
gender    0.006613642 0.08968717 -0.1673214 0.1817552
migration 0.006613642 0.08968717 -0.1673214 0.1817552

, , zi_persistence_dev1:anonymity_dev1

             Estimate  Est.Error       Q2.5     Q97.5
climate   0.008806348 0.08889244 -0.1644127 0.1844408
gender    0.008806348 0.08889244 -0.1644127 0.1844408
migration 0.008806348 0.08889244 -0.1644127 0.1844408


$`topic:group`
, , Intercept

                   Estimate Est.Error     Q2.5    Q97.5
climate_Madrid1    1.759102 0.2462320 1.239356 2.240928
climate_Madrid2    2.457696 0.2344562 1.953004 2.931168
climate_Madrid3    2.272688 0.2356103 1.761028 2.755240
climate_Madrid4    2.792250 0.2333683 2.290435 3.268066
climate_Oslo1      1.984540 0.2399568 1.473012 2.458905
climate_Oslo2      2.675209 0.2333150 2.186113 3.147148
climate_Oslo3      2.351983 0.2385897 1.854177 2.831614
climate_Oslo4      2.763715 0.2329267 2.260392 3.241184
climate_Prag1      2.238111 0.2407547 1.745460 2.709367
climate_Prag2      2.285453 0.2405552 1.779591 2.751184
climate_Prag3      2.799090 0.2332121 2.307405 3.258438
climate_Prag4      2.562857 0.2378152 2.068329 3.028966
climate_Rio1       2.393871 0.2370286 1.879201 2.869942
climate_Rio2       2.526772 0.2375326 2.019074 2.988996
climate_Rio3       2.564091 0.2333346 2.056815 3.018157
climate_Rio4       2.342742 0.2342666 1.835705 2.821173
gender_Dubai1      1.908363 0.2450134 1.391224 2.399666
gender_Dubai2      2.309364 0.2340309 1.815465 2.780017
gender_Dubai3      2.126913 0.2390267 1.631461 2.595698
gender_Dubai4      2.483903 0.2324739 1.993163 2.956572
gender_London1     2.295351 0.2327935 1.816049 2.776759
gender_London2     2.530184 0.2319773 2.054167 3.003108
gender_London3     3.291772 0.2277683 2.821858 3.762740
gender_London4     1.927815 0.2363204 1.441430 2.422490
gender_Rom1        2.478452 0.2313663 1.988047 2.964336
gender_Rom2        2.968613 0.2261618 2.482351 3.449909
gender_Rom3        2.114065 0.2429499 1.594716 2.612889
gender_Rom4        1.870895 0.2389636 1.378278 2.365217
gender_Tokio1      1.874823 0.2383897 1.370670 2.357965
gender_Tokio2      2.406254 0.2307610 1.914750 2.877738
gender_Tokio3      3.084374 0.2281562 2.598814 3.559181
gender_Tokio4      2.796354 0.2287945 2.301066 3.281697
migration_Berlin1  2.236043 0.2427117 1.775399 2.772168
migration_Berlin2  1.975575 0.2425887 1.512872 2.530394
migration_Berlin3  2.421517 0.2411588 1.957223 2.961730
migration_Berlin4  3.316509 0.2327748 2.875410 3.844989
migration_Florenz1 1.816609 0.2453584 1.330089 2.371305
migration_Florenz2 2.454008 0.2411424 1.982611 3.000546
migration_Florenz3 2.180384 0.2435751 1.720796 2.729374
migration_Florenz4 2.301909 0.2408945 1.850372 2.853787
migration_Paris1   1.973939 0.2450204 1.516412 2.509617
migration_Paris2   2.061176 0.2442268 1.598793 2.592717
migration_Paris3   2.271019 0.2386424 1.820523 2.806410
migration_Paris4   2.413627 0.2392892 1.963469 2.945058
migration_Sydney1  2.833173 0.2368622 2.389574 3.357629
migration_Sydney2  2.217617 0.2445668 1.761643 2.750194
migration_Sydney3  2.428997 0.2436732 1.970346 2.969473
migration_Sydney4  2.524944 0.2397629 2.073708 3.062403

, , zi_Intercept

                    Estimate Est.Error      Q2.5       Q97.5
climate_Madrid1    -1.370307 0.3721058 -2.106880 -0.61528568
climate_Madrid2    -1.421178 0.3827236 -2.191532 -0.66295948
climate_Madrid3    -1.314318 0.3680023 -2.006572 -0.55018676
climate_Madrid4    -1.314233 0.3615676 -1.988889 -0.55378347
climate_Oslo1      -1.268489 0.3648985 -1.928773 -0.49054187
climate_Oslo2      -1.271029 0.3629712 -1.919742 -0.50455122
climate_Oslo3      -1.322019 0.3658883 -2.020575 -0.55728205
climate_Oslo4      -1.116809 0.3840140 -1.724377 -0.24611549
climate_Prag1      -1.308407 0.3614901 -1.990049 -0.55732678
climate_Prag2      -1.205797 0.3640367 -1.827430 -0.39745096
climate_Prag3      -1.306936 0.3632909 -1.982392 -0.53620838
climate_Prag4      -1.307129 0.3618862 -1.977822 -0.56085374
climate_Rio1       -1.253380 0.3629223 -1.907779 -0.47604325
climate_Rio2       -1.149052 0.3774966 -1.772990 -0.29298313
climate_Rio3       -1.473682 0.4021293 -2.313559 -0.73098405
climate_Rio4       -1.307096 0.3653043 -1.987188 -0.54806953
gender_Dubai1      -1.098982 0.3848154 -1.695485 -0.21730384
gender_Dubai2      -1.355206 0.3690983 -2.070361 -0.61014911
gender_Dubai3      -1.356229 0.3746776 -2.082277 -0.58732582
gender_Dubai4      -1.356196 0.3705898 -2.089175 -0.60571297
gender_London1     -1.412310 0.3804841 -2.182613 -0.68963756
gender_London2     -1.303418 0.3665494 -1.993471 -0.54567104
gender_London3     -1.305189 0.3662344 -1.990580 -0.53238568
gender_London4     -1.361289 0.3683366 -2.084745 -0.60160100
gender_Rom1        -1.417191 0.3798392 -2.191277 -0.68170207
gender_Rom2        -1.207867 0.3660861 -1.834885 -0.38533330
gender_Rom3        -1.006564 0.4178990 -1.622727 -0.05748882
gender_Rom4        -1.259744 0.3649339 -1.908092 -0.47301603
gender_Tokio1      -1.422424 0.3831665 -2.210404 -0.68507654
gender_Tokio2      -1.316292 0.3642416 -2.005070 -0.55994924
gender_Tokio3      -1.211316 0.3647677 -1.831489 -0.40453563
gender_Tokio4      -1.368475 0.3721193 -2.108016 -0.61569715
migration_Berlin1  -1.362900 0.3711355 -2.088316 -0.60658403
migration_Berlin2  -1.365905 0.3711943 -2.094868 -0.63066407
migration_Berlin3  -1.105907 0.3877544 -1.705316 -0.20242532
migration_Berlin4  -1.412864 0.3808719 -2.187128 -0.67527560
migration_Florenz1 -1.373810 0.3775674 -2.153421 -0.61702415
migration_Florenz2 -1.262593 0.3630209 -1.923779 -0.48091321
migration_Florenz3 -1.312413 0.3633158 -1.984073 -0.55945309
migration_Florenz4 -1.316648 0.3661400 -2.010891 -0.55513195
migration_Paris1   -1.255209 0.3621782 -1.901247 -0.48386193
migration_Paris2   -1.147654 0.3765091 -1.753525 -0.29559249
migration_Paris3   -1.355056 0.3703160 -2.089092 -0.58375085
migration_Paris4   -1.304116 0.3643284 -1.981056 -0.54851894
migration_Sydney1  -1.410283 0.3816626 -2.182991 -0.66806776
migration_Sydney2  -1.145921 0.3780903 -1.764725 -0.28776382
migration_Sydney3  -1.250684 0.3656062 -1.900912 -0.44764680
migration_Sydney4  -1.409058 0.3825074 -2.186340 -0.64406722

, , persistence_dev1

                       Estimate  Est.Error       Q2.5     Q97.5
climate_Madrid1    -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Madrid2    -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Madrid3    -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Madrid4    -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Oslo1      -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Oslo2      -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Oslo3      -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Oslo4      -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Prag1      -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Prag2      -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Prag3      -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Prag4      -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Rio1       -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Rio2       -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Rio3       -0.005103491 0.05928069 -0.1215837 0.1109697
climate_Rio4       -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Dubai1      -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Dubai2      -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Dubai3      -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Dubai4      -0.005103491 0.05928069 -0.1215837 0.1109697
gender_London1     -0.005103491 0.05928069 -0.1215837 0.1109697
gender_London2     -0.005103491 0.05928069 -0.1215837 0.1109697
gender_London3     -0.005103491 0.05928069 -0.1215837 0.1109697
gender_London4     -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Rom1        -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Rom2        -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Rom3        -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Rom4        -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Tokio1      -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Tokio2      -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Tokio3      -0.005103491 0.05928069 -0.1215837 0.1109697
gender_Tokio4      -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Berlin1  -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Berlin2  -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Berlin3  -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Berlin4  -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Florenz1 -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Florenz2 -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Florenz3 -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Florenz4 -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Paris1   -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Paris2   -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Paris3   -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Paris4   -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Sydney1  -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Sydney2  -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Sydney3  -0.005103491 0.05928069 -0.1215837 0.1109697
migration_Sydney4  -0.005103491 0.05928069 -0.1215837 0.1109697

, , anonymity_dev1

                       Estimate  Est.Error       Q2.5     Q97.5
climate_Madrid1    0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Madrid2    0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Madrid3    0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Madrid4    0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Oslo1      0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Oslo2      0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Oslo3      0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Oslo4      0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Prag1      0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Prag2      0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Prag3      0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Prag4      0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Rio1       0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Rio2       0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Rio3       0.0009798597 0.05871981 -0.1159714 0.1169367
climate_Rio4       0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Dubai1      0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Dubai2      0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Dubai3      0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Dubai4      0.0009798597 0.05871981 -0.1159714 0.1169367
gender_London1     0.0009798597 0.05871981 -0.1159714 0.1169367
gender_London2     0.0009798597 0.05871981 -0.1159714 0.1169367
gender_London3     0.0009798597 0.05871981 -0.1159714 0.1169367
gender_London4     0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Rom1        0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Rom2        0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Rom3        0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Rom4        0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Tokio1      0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Tokio2      0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Tokio3      0.0009798597 0.05871981 -0.1159714 0.1169367
gender_Tokio4      0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Berlin1  0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Berlin2  0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Berlin3  0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Berlin4  0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Florenz1 0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Florenz2 0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Florenz3 0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Florenz4 0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Paris1   0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Paris2   0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Paris3   0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Paris4   0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Sydney1  0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Sydney2  0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Sydney3  0.0009798597 0.05871981 -0.1159714 0.1169367
migration_Sydney4  0.0009798597 0.05871981 -0.1159714 0.1169367

, , persistence_dev1:anonymity_dev1

                     Estimate  Est.Error        Q2.5     Q97.5
climate_Madrid1    0.02947201 0.05854357 -0.08580015 0.1467483
climate_Madrid2    0.02947201 0.05854357 -0.08580015 0.1467483
climate_Madrid3    0.02947201 0.05854357 -0.08580015 0.1467483
climate_Madrid4    0.02947201 0.05854357 -0.08580015 0.1467483
climate_Oslo1      0.02947201 0.05854357 -0.08580015 0.1467483
climate_Oslo2      0.02947201 0.05854357 -0.08580015 0.1467483
climate_Oslo3      0.02947201 0.05854357 -0.08580015 0.1467483
climate_Oslo4      0.02947201 0.05854357 -0.08580015 0.1467483
climate_Prag1      0.02947201 0.05854357 -0.08580015 0.1467483
climate_Prag2      0.02947201 0.05854357 -0.08580015 0.1467483
climate_Prag3      0.02947201 0.05854357 -0.08580015 0.1467483
climate_Prag4      0.02947201 0.05854357 -0.08580015 0.1467483
climate_Rio1       0.02947201 0.05854357 -0.08580015 0.1467483
climate_Rio2       0.02947201 0.05854357 -0.08580015 0.1467483
climate_Rio3       0.02947201 0.05854357 -0.08580015 0.1467483
climate_Rio4       0.02947201 0.05854357 -0.08580015 0.1467483
gender_Dubai1      0.02947201 0.05854357 -0.08580015 0.1467483
gender_Dubai2      0.02947201 0.05854357 -0.08580015 0.1467483
gender_Dubai3      0.02947201 0.05854357 -0.08580015 0.1467483
gender_Dubai4      0.02947201 0.05854357 -0.08580015 0.1467483
gender_London1     0.02947201 0.05854357 -0.08580015 0.1467483
gender_London2     0.02947201 0.05854357 -0.08580015 0.1467483
gender_London3     0.02947201 0.05854357 -0.08580015 0.1467483
gender_London4     0.02947201 0.05854357 -0.08580015 0.1467483
gender_Rom1        0.02947201 0.05854357 -0.08580015 0.1467483
gender_Rom2        0.02947201 0.05854357 -0.08580015 0.1467483
gender_Rom3        0.02947201 0.05854357 -0.08580015 0.1467483
gender_Rom4        0.02947201 0.05854357 -0.08580015 0.1467483
gender_Tokio1      0.02947201 0.05854357 -0.08580015 0.1467483
gender_Tokio2      0.02947201 0.05854357 -0.08580015 0.1467483
gender_Tokio3      0.02947201 0.05854357 -0.08580015 0.1467483
gender_Tokio4      0.02947201 0.05854357 -0.08580015 0.1467483
migration_Berlin1  0.02947201 0.05854357 -0.08580015 0.1467483
migration_Berlin2  0.02947201 0.05854357 -0.08580015 0.1467483
migration_Berlin3  0.02947201 0.05854357 -0.08580015 0.1467483
migration_Berlin4  0.02947201 0.05854357 -0.08580015 0.1467483
migration_Florenz1 0.02947201 0.05854357 -0.08580015 0.1467483
migration_Florenz2 0.02947201 0.05854357 -0.08580015 0.1467483
migration_Florenz3 0.02947201 0.05854357 -0.08580015 0.1467483
migration_Florenz4 0.02947201 0.05854357 -0.08580015 0.1467483
migration_Paris1   0.02947201 0.05854357 -0.08580015 0.1467483
migration_Paris2   0.02947201 0.05854357 -0.08580015 0.1467483
migration_Paris3   0.02947201 0.05854357 -0.08580015 0.1467483
migration_Paris4   0.02947201 0.05854357 -0.08580015 0.1467483
migration_Sydney1  0.02947201 0.05854357 -0.08580015 0.1467483
migration_Sydney2  0.02947201 0.05854357 -0.08580015 0.1467483
migration_Sydney3  0.02947201 0.05854357 -0.08580015 0.1467483
migration_Sydney4  0.02947201 0.05854357 -0.08580015 0.1467483

, , zi_persistence_dev1

                     Estimate  Est.Error      Q2.5     Q97.5
climate_Madrid1    0.02959441 0.08905165 -0.148826 0.2034417
climate_Madrid2    0.02959441 0.08905165 -0.148826 0.2034417
climate_Madrid3    0.02959441 0.08905165 -0.148826 0.2034417
climate_Madrid4    0.02959441 0.08905165 -0.148826 0.2034417
climate_Oslo1      0.02959441 0.08905165 -0.148826 0.2034417
climate_Oslo2      0.02959441 0.08905165 -0.148826 0.2034417
climate_Oslo3      0.02959441 0.08905165 -0.148826 0.2034417
climate_Oslo4      0.02959441 0.08905165 -0.148826 0.2034417
climate_Prag1      0.02959441 0.08905165 -0.148826 0.2034417
climate_Prag2      0.02959441 0.08905165 -0.148826 0.2034417
climate_Prag3      0.02959441 0.08905165 -0.148826 0.2034417
climate_Prag4      0.02959441 0.08905165 -0.148826 0.2034417
climate_Rio1       0.02959441 0.08905165 -0.148826 0.2034417
climate_Rio2       0.02959441 0.08905165 -0.148826 0.2034417
climate_Rio3       0.02959441 0.08905165 -0.148826 0.2034417
climate_Rio4       0.02959441 0.08905165 -0.148826 0.2034417
gender_Dubai1      0.02959441 0.08905165 -0.148826 0.2034417
gender_Dubai2      0.02959441 0.08905165 -0.148826 0.2034417
gender_Dubai3      0.02959441 0.08905165 -0.148826 0.2034417
gender_Dubai4      0.02959441 0.08905165 -0.148826 0.2034417
gender_London1     0.02959441 0.08905165 -0.148826 0.2034417
gender_London2     0.02959441 0.08905165 -0.148826 0.2034417
gender_London3     0.02959441 0.08905165 -0.148826 0.2034417
gender_London4     0.02959441 0.08905165 -0.148826 0.2034417
gender_Rom1        0.02959441 0.08905165 -0.148826 0.2034417
gender_Rom2        0.02959441 0.08905165 -0.148826 0.2034417
gender_Rom3        0.02959441 0.08905165 -0.148826 0.2034417
gender_Rom4        0.02959441 0.08905165 -0.148826 0.2034417
gender_Tokio1      0.02959441 0.08905165 -0.148826 0.2034417
gender_Tokio2      0.02959441 0.08905165 -0.148826 0.2034417
gender_Tokio3      0.02959441 0.08905165 -0.148826 0.2034417
gender_Tokio4      0.02959441 0.08905165 -0.148826 0.2034417
migration_Berlin1  0.02959441 0.08905165 -0.148826 0.2034417
migration_Berlin2  0.02959441 0.08905165 -0.148826 0.2034417
migration_Berlin3  0.02959441 0.08905165 -0.148826 0.2034417
migration_Berlin4  0.02959441 0.08905165 -0.148826 0.2034417
migration_Florenz1 0.02959441 0.08905165 -0.148826 0.2034417
migration_Florenz2 0.02959441 0.08905165 -0.148826 0.2034417
migration_Florenz3 0.02959441 0.08905165 -0.148826 0.2034417
migration_Florenz4 0.02959441 0.08905165 -0.148826 0.2034417
migration_Paris1   0.02959441 0.08905165 -0.148826 0.2034417
migration_Paris2   0.02959441 0.08905165 -0.148826 0.2034417
migration_Paris3   0.02959441 0.08905165 -0.148826 0.2034417
migration_Paris4   0.02959441 0.08905165 -0.148826 0.2034417
migration_Sydney1  0.02959441 0.08905165 -0.148826 0.2034417
migration_Sydney2  0.02959441 0.08905165 -0.148826 0.2034417
migration_Sydney3  0.02959441 0.08905165 -0.148826 0.2034417
migration_Sydney4  0.02959441 0.08905165 -0.148826 0.2034417

, , zi_anonymity_dev1

                      Estimate  Est.Error       Q2.5     Q97.5
climate_Madrid1    0.006613642 0.08968717 -0.1673214 0.1817552
climate_Madrid2    0.006613642 0.08968717 -0.1673214 0.1817552
climate_Madrid3    0.006613642 0.08968717 -0.1673214 0.1817552
climate_Madrid4    0.006613642 0.08968717 -0.1673214 0.1817552
climate_Oslo1      0.006613642 0.08968717 -0.1673214 0.1817552
climate_Oslo2      0.006613642 0.08968717 -0.1673214 0.1817552
climate_Oslo3      0.006613642 0.08968717 -0.1673214 0.1817552
climate_Oslo4      0.006613642 0.08968717 -0.1673214 0.1817552
climate_Prag1      0.006613642 0.08968717 -0.1673214 0.1817552
climate_Prag2      0.006613642 0.08968717 -0.1673214 0.1817552
climate_Prag3      0.006613642 0.08968717 -0.1673214 0.1817552
climate_Prag4      0.006613642 0.08968717 -0.1673214 0.1817552
climate_Rio1       0.006613642 0.08968717 -0.1673214 0.1817552
climate_Rio2       0.006613642 0.08968717 -0.1673214 0.1817552
climate_Rio3       0.006613642 0.08968717 -0.1673214 0.1817552
climate_Rio4       0.006613642 0.08968717 -0.1673214 0.1817552
gender_Dubai1      0.006613642 0.08968717 -0.1673214 0.1817552
gender_Dubai2      0.006613642 0.08968717 -0.1673214 0.1817552
gender_Dubai3      0.006613642 0.08968717 -0.1673214 0.1817552
gender_Dubai4      0.006613642 0.08968717 -0.1673214 0.1817552
gender_London1     0.006613642 0.08968717 -0.1673214 0.1817552
gender_London2     0.006613642 0.08968717 -0.1673214 0.1817552
gender_London3     0.006613642 0.08968717 -0.1673214 0.1817552
gender_London4     0.006613642 0.08968717 -0.1673214 0.1817552
gender_Rom1        0.006613642 0.08968717 -0.1673214 0.1817552
gender_Rom2        0.006613642 0.08968717 -0.1673214 0.1817552
gender_Rom3        0.006613642 0.08968717 -0.1673214 0.1817552
gender_Rom4        0.006613642 0.08968717 -0.1673214 0.1817552
gender_Tokio1      0.006613642 0.08968717 -0.1673214 0.1817552
gender_Tokio2      0.006613642 0.08968717 -0.1673214 0.1817552
gender_Tokio3      0.006613642 0.08968717 -0.1673214 0.1817552
gender_Tokio4      0.006613642 0.08968717 -0.1673214 0.1817552
migration_Berlin1  0.006613642 0.08968717 -0.1673214 0.1817552
migration_Berlin2  0.006613642 0.08968717 -0.1673214 0.1817552
migration_Berlin3  0.006613642 0.08968717 -0.1673214 0.1817552
migration_Berlin4  0.006613642 0.08968717 -0.1673214 0.1817552
migration_Florenz1 0.006613642 0.08968717 -0.1673214 0.1817552
migration_Florenz2 0.006613642 0.08968717 -0.1673214 0.1817552
migration_Florenz3 0.006613642 0.08968717 -0.1673214 0.1817552
migration_Florenz4 0.006613642 0.08968717 -0.1673214 0.1817552
migration_Paris1   0.006613642 0.08968717 -0.1673214 0.1817552
migration_Paris2   0.006613642 0.08968717 -0.1673214 0.1817552
migration_Paris3   0.006613642 0.08968717 -0.1673214 0.1817552
migration_Paris4   0.006613642 0.08968717 -0.1673214 0.1817552
migration_Sydney1  0.006613642 0.08968717 -0.1673214 0.1817552
migration_Sydney2  0.006613642 0.08968717 -0.1673214 0.1817552
migration_Sydney3  0.006613642 0.08968717 -0.1673214 0.1817552
migration_Sydney4  0.006613642 0.08968717 -0.1673214 0.1817552

, , zi_persistence_dev1:anonymity_dev1

                      Estimate  Est.Error       Q2.5     Q97.5
climate_Madrid1    0.008806348 0.08889244 -0.1644127 0.1844408
climate_Madrid2    0.008806348 0.08889244 -0.1644127 0.1844408
climate_Madrid3    0.008806348 0.08889244 -0.1644127 0.1844408
climate_Madrid4    0.008806348 0.08889244 -0.1644127 0.1844408
climate_Oslo1      0.008806348 0.08889244 -0.1644127 0.1844408
climate_Oslo2      0.008806348 0.08889244 -0.1644127 0.1844408
climate_Oslo3      0.008806348 0.08889244 -0.1644127 0.1844408
climate_Oslo4      0.008806348 0.08889244 -0.1644127 0.1844408
climate_Prag1      0.008806348 0.08889244 -0.1644127 0.1844408
climate_Prag2      0.008806348 0.08889244 -0.1644127 0.1844408
climate_Prag3      0.008806348 0.08889244 -0.1644127 0.1844408
climate_Prag4      0.008806348 0.08889244 -0.1644127 0.1844408
climate_Rio1       0.008806348 0.08889244 -0.1644127 0.1844408
climate_Rio2       0.008806348 0.08889244 -0.1644127 0.1844408
climate_Rio3       0.008806348 0.08889244 -0.1644127 0.1844408
climate_Rio4       0.008806348 0.08889244 -0.1644127 0.1844408
gender_Dubai1      0.008806348 0.08889244 -0.1644127 0.1844408
gender_Dubai2      0.008806348 0.08889244 -0.1644127 0.1844408
gender_Dubai3      0.008806348 0.08889244 -0.1644127 0.1844408
gender_Dubai4      0.008806348 0.08889244 -0.1644127 0.1844408
gender_London1     0.008806348 0.08889244 -0.1644127 0.1844408
gender_London2     0.008806348 0.08889244 -0.1644127 0.1844408
gender_London3     0.008806348 0.08889244 -0.1644127 0.1844408
gender_London4     0.008806348 0.08889244 -0.1644127 0.1844408
gender_Rom1        0.008806348 0.08889244 -0.1644127 0.1844408
gender_Rom2        0.008806348 0.08889244 -0.1644127 0.1844408
gender_Rom3        0.008806348 0.08889244 -0.1644127 0.1844408
gender_Rom4        0.008806348 0.08889244 -0.1644127 0.1844408
gender_Tokio1      0.008806348 0.08889244 -0.1644127 0.1844408
gender_Tokio2      0.008806348 0.08889244 -0.1644127 0.1844408
gender_Tokio3      0.008806348 0.08889244 -0.1644127 0.1844408
gender_Tokio4      0.008806348 0.08889244 -0.1644127 0.1844408
migration_Berlin1  0.008806348 0.08889244 -0.1644127 0.1844408
migration_Berlin2  0.008806348 0.08889244 -0.1644127 0.1844408
migration_Berlin3  0.008806348 0.08889244 -0.1644127 0.1844408
migration_Berlin4  0.008806348 0.08889244 -0.1644127 0.1844408
migration_Florenz1 0.008806348 0.08889244 -0.1644127 0.1844408
migration_Florenz2 0.008806348 0.08889244 -0.1644127 0.1844408
migration_Florenz3 0.008806348 0.08889244 -0.1644127 0.1844408
migration_Florenz4 0.008806348 0.08889244 -0.1644127 0.1844408
migration_Paris1   0.008806348 0.08889244 -0.1644127 0.1844408
migration_Paris2   0.008806348 0.08889244 -0.1644127 0.1844408
migration_Paris3   0.008806348 0.08889244 -0.1644127 0.1844408
migration_Paris4   0.008806348 0.08889244 -0.1644127 0.1844408
migration_Sydney1  0.008806348 0.08889244 -0.1644127 0.1844408
migration_Sydney2  0.008806348 0.08889244 -0.1644127 0.1844408
migration_Sydney3  0.008806348 0.08889244 -0.1644127 0.1844408
migration_Sydney4  0.008806348 0.08889244 -0.1644127 0.1844408

Let’s look at exponentiated results for interpretation.

broom.mixed::tidy(fit_hrdl_1) |> 
  mutate(
    estimate_exp = exp(estimate),
    conf_low_exp = exp(conf.low),
    conf_high_exp = exp(conf.high),
    estimate_zi = exp(estimate) / (1 + exp(estimate)),
    conf_low_zi = exp(conf.low) / (1 + exp(conf.low)),
    conf_high_zi = exp(conf.high) / (1 + exp(conf.high))
    ) %>% 
  mutate(
    across(
      c(
        estimate
        , conf.low
        , conf.high
        , estimate_exp
        , conf_low_exp
        , conf_high_exp
        , estimate_zi 
        , conf_low_zi 
        , conf_high_zi
        ),
      ~ round(.x, 2)
    )
  ) %>% 
  select(
    -effect, -group, -std.error  # -component, 
  ) %>% 
  rmarkdown::paged_table()
Warning in tidy.brmsfit(fit_hrdl_1): some parameter names contain underscores:
term naming may be unreliable!

Let’s visualize results.

plot(
  conditional_effects(
    fit_hrdl_1
    ), 
  ask = FALSE
  )

Similar picture. Effects appear even smaller.

Let’s look at posteriors

p_4 <- 
  pp_check(fit_hrdl_1) + 
  labs(title = "Zero-inflated poisson")
Using 10 posterior draws for ppc type 'dens_overlay' by default.
p_4

Exploratory Analyses

Frequentist

Look at results from a frequentist perspective.

Fixed effects

Estimate nested model.

fit_fe_1_frq <- 
  lmer(
    op_expr ~ 
      1 + 
      (1 | topic/group) + 
      persistence_dev * anonymity_dev
    , data = d
    )

summary(fit_fe_1_frq)
Linear mixed model fit by REML ['lmerMod']
Formula: op_expr ~ 1 + (1 | topic/group) + persistence_dev * anonymity_dev
   Data: d

REML criterion at convergence: 7606.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4108 -0.5454 -0.2147  0.2659 13.5411 

Random effects:
 Groups      Name        Variance Std.Dev.
 group:topic (Intercept)  10.94    3.307  
 topic       (Intercept)   0.00    0.000  
 Residual                155.90   12.486  
Number of obs: 960, groups:  group:topic, 48; topic, 3

Fixed effects:
                                Estimate Std. Error t value
(Intercept)                       9.2729     0.6247  14.844
persistence_dev1                 -0.0125     0.6247  -0.020
anonymity_dev1                   -0.2417     0.6247  -0.387
persistence_dev1:anonymity_dev1   0.2521     0.6247   0.404

Correlation of Fixed Effects:
            (Intr) prss_1 anny_1
prsstnc_dv1 0.000               
annymty_dv1 0.000  0.000        
prsstn_1:_1 0.000  0.000  0.000 
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Quite weird that topic doesn’t get any variance at all. Perhaps due to small cluster size? With Bayesian estimation, it worked alright.

Estimate without nesting.

fit_fe_2_frq <- 
  lmer(
    op_expr ~ 
      1 + 
      (1 | topic) +
      persistence_dev * anonymity_dev
    , data = d
    )

summary(fit_fe_2_frq)
Linear mixed model fit by REML ['lmerMod']
Formula: op_expr ~ 1 + (1 | topic) + persistence_dev * anonymity_dev
   Data: d

REML criterion at convergence: 7627

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.7806 -0.6106 -0.2378  0.2292 13.7610 

Random effects:
 Groups   Name        Variance Std.Dev.
 topic    (Intercept)   0.3331  0.5771 
 Residual             165.7407 12.8740 
Number of obs: 960, groups:  topic, 3

Fixed effects:
                                Estimate Std. Error t value
(Intercept)                       9.2729     0.5326  17.410
persistence_dev1                 -0.0125     0.4155  -0.030
anonymity_dev1                   -0.2417     0.4155  -0.582
persistence_dev1:anonymity_dev1   0.2521     0.4155   0.607

Correlation of Fixed Effects:
            (Intr) prss_1 anny_1
prsstnc_dv1 0.000               
annymty_dv1 0.000  0.000        
prsstn_1:_1 0.000  0.000  0.000 

Estimate without hierarchical structure.

fit_fe_3_frq <- 
  lm(
    op_expr ~ 
      1 + 
      persistence_dev * anonymity_dev + topic
    , data = d
    )

summary(fit_fe_3_frq)

Call:
lm(formula = op_expr ~ 1 + persistence_dev * anonymity_dev + 
    topic, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.469  -7.744  -3.150   3.040 177.798 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       9.6312     0.7197  13.383   <2e-16 ***
persistence_dev1                 -0.0125     0.4155  -0.030    0.976    
anonymity_dev1                   -0.2417     0.4155  -0.582    0.561    
topicgender                       0.3312     1.0178   0.325    0.745    
topicmigration                   -1.4062     1.0178  -1.382    0.167    
persistence_dev1:anonymity_dev1   0.2521     0.4155   0.607    0.544    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.87 on 954 degrees of freedom
Multiple R-squared:  0.004169,  Adjusted R-squared:  -0.001051 
F-statistic: 0.7987 on 5 and 954 DF,  p-value: 0.5507

Also here, no significant effects.

Random effects

Now do the same with random effects.

fit_re_1_frq <- 
  lmer(
    op_expr ~ 
      1 + 
      persistence_dev * anonymity_dev + 
      (1 + persistence_dev * anonymity_dev | topic) +
      (1 | topic:group),
    , data = d
    )

summary(fit_re_1_frq)
Linear mixed model fit by REML ['lmerMod']
Formula: 
op_expr ~ 1 + persistence_dev * anonymity_dev + (1 + persistence_dev *  
    anonymity_dev | topic) + (1 | topic:group)
   Data: d

REML criterion at convergence: 7603.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4310 -0.5485 -0.2174  0.2635 13.5713 

Random effects:
 Groups      Name                            Variance  Std.Dev. Corr       
 topic:group (Intercept)                     8.751e+00  2.95816            
 topic       (Intercept)                     5.858e-01  0.76535            
             persistence_dev1                1.752e-02  0.13238 -1.00      
             anonymity_dev1                  8.313e-03  0.09117 -1.00  1.00
             persistence_dev1:anonymity_dev1 2.396e+00  1.54785  1.00 -1.00
 Residual                                    1.559e+02 12.48581            
      
      
      
      
      
 -1.00
      
Number of obs: 960, groups:  topic:group, 48; topic, 3

Fixed effects:
                                Estimate Std. Error t value
(Intercept)                       9.2729     0.7348  12.619
persistence_dev1                 -0.0125     0.5921  -0.021
anonymity_dev1                   -0.2417     0.5895  -0.410
persistence_dev1:anonymity_dev1   0.2521     1.0693   0.236

Correlation of Fixed Effects:
            (Intr) prss_1 anny_1
prsstnc_dv1 -0.078              
annymty_dv1 -0.054  0.012       
prsstn_1:_1  0.503 -0.108 -0.075
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
fit_re_2_frq <- 
  lmer(
    op_expr ~ 
      1 + 
      persistence_dev * anonymity_dev + 
      (1 + persistence_dev * anonymity_dev | topic),
    , data = d
    )

summary(fit_re_2_frq)
Linear mixed model fit by REML ['lmerMod']
Formula: 
op_expr ~ 1 + persistence_dev * anonymity_dev + (1 + persistence_dev *  
    anonymity_dev | topic)
   Data: d

REML criterion at convergence: 7617.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.9574 -0.5665 -0.2408  0.2263 13.7512 

Random effects:
 Groups   Name                            Variance  Std.Dev. Corr             
 topic    (Intercept)                     6.878e-01  0.82936                  
          persistence_dev1                2.061e-02  0.14355 -1.00            
          anonymity_dev1                  9.755e-03  0.09876 -1.00  1.00      
          persistence_dev1:anonymity_dev1 2.814e+00  1.67745  1.00 -1.00 -1.00
 Residual                                 1.636e+02 12.79063                  
Number of obs: 960, groups:  topic, 3

Fixed effects:
                                Estimate Std. Error t value
(Intercept)                       9.2729     0.6322  14.667
persistence_dev1                 -0.0125     0.4210  -0.030
anonymity_dev1                   -0.2417     0.4167  -0.580
persistence_dev1:anonymity_dev1   0.2521     1.0528   0.239

Correlation of Fixed Effects:
            (Intr) prss_1 anny_1
prsstnc_dv1 -0.149              
annymty_dv1 -0.104  0.027       
prsstn_1:_1  0.697 -0.181 -0.126
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Let’s inspect coefficients.

coefficients(fit_re_2_frq)
$topic
          (Intercept) persistence_dev1 anonymity_dev1
climate      9.331303      -0.02260573     -0.2486197
gender      10.017024      -0.14129197     -0.3302795
migration    8.470422       0.12639770     -0.1461008
          persistence_dev1:anonymity_dev1
climate                         0.3701762
gender                          1.7571129
migration                      -1.3710392

attr(,"class")
[1] "coef.mer"

Gender

Let’s see if effects differ across genders

fit_fe_gen_1 <- 
  brm(
    op_expr ~ 
      1 + persistence_dev * anonymity_dev +
      (1 | topic/group)
    , data = d
    , chains = 4
    , cores = 4
    , iter = 8000
    , warmup = 2000
    , family = zero_inflated_poisson()
    , control = list(
      adapt_delta = .95
      , max_treedepth = 12
      )
    )
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
using SDK: ‘MacOSX14.4.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Warning: There were 315 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems

Shows some problems, but not all that many.

Let’s inspect model.

plot(fit_fe_gen_1, ask = FALSE)

Traceplots look alright.

Let’s look at results.

summary(fit_fe_gen_1)
Warning: There were 315 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
 Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev + (1 | topic/group) 
   Data: d (Number of observations: 960) 
  Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
         total post-warmup draws = 24000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.27      0.33     0.01     1.21 1.00     2388     3643

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.40      0.05     0.32     0.50 1.00     3576     5264

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           2.39      0.20     1.95     2.82 1.00
persistence_dev1                   -0.00      0.06    -0.12     0.11 1.00
anonymity_dev1                      0.00      0.06    -0.12     0.12 1.00
persistence_dev1:anonymity_dev1     0.03      0.06    -0.09     0.15 1.00
                                Bulk_ESS Tail_ESS
Intercept                           4302     2993
persistence_dev1                    4061     6358
anonymity_dev1                      4270     7001
persistence_dev1:anonymity_dev1     3856     5126

Further Distributional Parameters:
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi     0.21      0.01     0.19     0.24 1.00    16100    13653

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Benefits

Let’s see if benefits differ across experimental groups.

We first look at the experimental group’s descriptives

d %>% 
  group_by(persistence) %>% 
  summarize(benefits_m = mean(benefits, na.rm = TRUE))
# A tibble: 2 × 2
  persistence benefits_m
  <chr>            <dbl>
1 high              3.12
2 low               3.23

Looking at persistence, we see people with lower persistence reporting slightly higher benefits.

d %>% 
  group_by(anonymity) %>% 
  summarize(benefits_m = mean(benefits, na.rm = TRUE))
# A tibble: 2 × 2
  anonymity benefits_m
  <chr>          <dbl>
1 high            3.14
2 low             3.20

Looking at anonymity, we see people with low anonymity reporting marginally higher benefits.

d %>% 
  group_by(persistence, anonymity) %>% 
  summarize(benefits_m = mean(benefits, na.rm = T))
`summarise()` has grouped output by 'persistence'. You can override using the
`.groups` argument.
# A tibble: 4 × 3
# Groups:   persistence [2]
  persistence anonymity benefits_m
  <chr>       <chr>          <dbl>
1 high        high            3.07
2 high        low             3.18
3 low         high            3.22
4 low         low             3.23

Looking at both groups combined, we see that low anonymity and low persistence yielded highest benefits.

fit_fe_ben_1 <- 
  brm(
    benefits ~ 
      1 + persistence_dev * anonymity_dev +
      (1 | topic/group)
    , data = d
    , chains = 4
    , cores = 4
    , iter = 8000
    , warmup = 2000
    , control = list(
      adapt_delta = .95
      , max_treedepth = 12
      )
    )
Warning: Rows containing NAs were excluded from the model.
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
using SDK: ‘MacOSX14.4.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Warning: There were 146 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems

Shows some problems, but not all that many.

Let’s inspect model.

plot(fit_fe_ben_1, ask = FALSE)

Traceplots look alright.

Let’s look at results.

summary(fit_fe_ben_1)
Warning: There were 146 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: benefits ~ 1 + persistence_dev * anonymity_dev + (1 | topic/group) 
   Data: d (Number of observations: 705) 
  Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
         total post-warmup draws = 24000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.14      0.22     0.00     0.78 1.00     3274     2146

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.07      0.05     0.00     0.17 1.00     6912     9608

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           3.18      0.11     2.93     3.44 1.00
persistence_dev1                   -0.05      0.03    -0.11     0.01 1.00
anonymity_dev1                     -0.03      0.03    -0.09     0.03 1.00
persistence_dev1:anonymity_dev1    -0.02      0.03    -0.08     0.04 1.00
                                Bulk_ESS Tail_ESS
Intercept                           5292     2048
persistence_dev1                   25336    17342
anonymity_dev1                     24256    16751
persistence_dev1:anonymity_dev1    24260    15633

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.74      0.02     0.70     0.78 1.00    29448    16416

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

No significant effect. But note that effect of persistence on perceived benefits only marginally not significant.

Risks

Let’s see if perceived differed across experimental groups.

We first look at the experimental group’s descriptives

d %>% 
  group_by(persistence) %>% 
  summarize(costs = mean(costs, na.rm = TRUE))
# A tibble: 2 × 2
  persistence costs
  <chr>       <dbl>
1 high         1.99
2 low          1.99

Looking at persistence, we see both groups report equal costs.

d %>% 
  group_by(anonymity) %>% 
  summarize(costs = mean(costs, na.rm = TRUE))
# A tibble: 2 × 2
  anonymity costs
  <chr>     <dbl>
1 high       1.89
2 low        2.09

Looking at anonymity, we see people with low anonymity report slightly higher costs.

d %>% 
  group_by(persistence, anonymity) %>% 
  summarize(costs = mean(costs, na.rm = TRUE))
`summarise()` has grouped output by 'persistence'. You can override using the
`.groups` argument.
# A tibble: 4 × 3
# Groups:   persistence [2]
  persistence anonymity costs
  <chr>       <chr>     <dbl>
1 high        high       1.90
2 high        low        2.07
3 low         high       1.87
4 low         low        2.11

Looking at both groups combined, we see that highest costs were reported by group with low anonymity and low persistence.

fit_fe_ris_1 <- 
  brm(
    costs ~ 
      1 + persistence_dev * anonymity_dev +
      (1 | topic/group)
    , data = d
    , chains = 4
    , cores = 4
    , iter = 20000
    , warmup = 2000
    , control = list(
      adapt_delta = .95
      , max_treedepth = 12
      )
    )
Warning: Rows containing NAs were excluded from the model.
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
using SDK: ‘MacOSX14.4.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Warning: There were 495 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems

Shows some problems, but not all that many.

Let’s inspect model.

plot(fit_fe_ris_1, ask = FALSE)

Traceplots look alright.

Let’s look at results.

summary(fit_fe_ris_1)
Warning: There were 495 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: costs ~ 1 + persistence_dev * anonymity_dev + (1 | topic/group) 
   Data: d (Number of observations: 705) 
  Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
         total post-warmup draws = 72000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.14      0.21     0.00     0.74 1.00    11201     5353

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.07      0.05     0.00     0.17 1.00    19946    24053

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           1.99      0.11     1.74     2.23 1.00
persistence_dev1                   -0.00      0.03    -0.07     0.06 1.00
anonymity_dev1                     -0.10      0.03    -0.17    -0.03 1.00
persistence_dev1:anonymity_dev1     0.02      0.03    -0.05     0.08 1.00
                                Bulk_ESS Tail_ESS
Intercept                          10504     5540
persistence_dev1                   62661    47237
anonymity_dev1                     65846    48554
persistence_dev1:anonymity_dev1    58363    30821

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.85      0.02     0.81     0.90 1.00    75075    46472

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

We find that anonymity does reduce costs.

Mediation

Let’s see if perceived benefits and risks were associated with increased opinion expressions.

fit_fe_med <- 
  brm(
    op_expr ~ 
      1 + persistence_dev * anonymity_dev + benefits + costs +
      (1 | topic/group)
    , data = d
    , chains = 4
    , cores = 4
    , iter = 8000
    , warmup = 2000
    , family = zero_inflated_poisson()
    , control = list(
      adapt_delta = .95
      , max_treedepth = 12
      )
    )
Warning: Rows containing NAs were excluded from the model.
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
using SDK: ‘MacOSX14.4.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Start sampling
Warning: There were 351 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems

Let’s look at results.

summary(fit_fe_med)
Warning: There were 351 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
 Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev + benefits + costs + (1 | topic/group) 
   Data: d (Number of observations: 705) 
  Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
         total post-warmup draws = 24000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.21      0.25     0.01     0.96 1.00     3142     7447

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.41      0.05     0.33     0.51 1.00     4402     7783

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           2.31      0.16     1.94     2.66 1.00
persistence_dev1                   -0.02      0.06    -0.14     0.10 1.00
anonymity_dev1                      0.01      0.06    -0.11     0.13 1.00
benefits                            0.11      0.02     0.08     0.14 1.00
costs                              -0.12      0.01    -0.14    -0.09 1.00
persistence_dev1:anonymity_dev1     0.03      0.06    -0.09     0.15 1.00
                                Bulk_ESS Tail_ESS
Intercept                           6722     5749
persistence_dev1                    3874     6750
anonymity_dev1                      3995     6829
benefits                           16490    11997
costs                              16979    15573
persistence_dev1:anonymity_dev1     3881     7316

Further Distributional Parameters:
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi     0.08      0.01     0.06     0.10 1.00    17823    13042

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

We find that increased perceived costs are associated with decreased opinion expressions. Increased benefits are associated with increased opinion expressions. Let’s check if overall effect is significant.

anon_risk_a_b <- fixef(fit_fe_ris_1)["anonymity_dev1", "Estimate"]
anon_risk_a_se <- fixef(fit_fe_ris_1)["anonymity_dev1", "Est.Error"]
anon_risk_a_dis <- rnorm(10000, anon_risk_a_b, anon_risk_a_se)

anon_risk_b_b <- fixef(fit_fe_med)["benefits", "Estimate"]
anon_risk_b_se <- fixef(fit_fe_med)["benefits", "Est.Error"]
anon_risk_b_dis <- rnorm(10000, anon_risk_b_b, anon_risk_b_se)

anon_risk_ab_dis <- anon_risk_a_dis * anon_risk_b_dis
anon_risk_ab_m <- median(anon_risk_ab_dis)
anon_risk_ab_ll <- quantile(anon_risk_ab_dis, .025)
anon_risk_ab_ul <- quantile(anon_risk_ab_dis, .975)

The effect is significant (b = -0.0107919, 95% MC CI [-0.0198564, -0.0034263]).

Save

save.image("data/image.RData")